

use book-basedata-replication, clear

*estimate model
quietly stcox med2 med2_t prevcris2 viol2 crisdur2 jointdem victory2 contig2 /* 
*/		prevcris2_t viol2_t crisdur2_t jointdem_t victory2_t contig2_t /* 
*/		if _t<3650,  strata(order5) cluster(dyadno) efron nohr


*calculate mean of each variable
foreach var of varlist med2 prevcris2 viol2 crisdur2 jointdem victory2 contig2 {
quietly sum `var' if e(sample)
local s_`var'=r(mean)
display `s_`var''
}

*calculate survival curves
scurve_tvc if _t<3650, gen(base_S_mean) at(med2 0 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)
scurve_tvc if _t<3650, gen(med_S_mean) at(med2 1 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)

*calculate difference in survival curves for each stratum
forvalues i=1/5 {
gen Sdiff`i'=med_S_mean`i'-base_S_mean`i'
}

*run bootstrap calculations 
preserve 
clear
*generate variables to store bootstrap results
forvalues i=1/5 {
gen S_diff`i'=.
}
gen _tscurve=.
save data.dta, replace
restore

*calculate 1000 replications
forvalues i=1/1000 {
preserve
*generate bootstrap sample
bsample, cluster(id)
*calculate survival curves for bootstrap sample
scurve_tvc if _t<3650, gen(S0_) at(med2 0 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)
scurve_tvc if _t<3650, gen(S1_) at(med2 1 prevcris2 `s_prevcris2' /* 
*/		viol2 `s_viol2' crisdur2 `s_crisdur2' jointdem `s_jointdem' /* 
*/		victory2 `s_victory2' contig2 `s_contig2') tvc(med2) texp(_t) /*
*/ 		strata(order5) ties(efron)

*calculate difference in survival curves in each stratum for bootstrap sample
forvalues j=1/5 {
gen S_diff`j'=S1_`j'-S0_`j'
}

*keep difference and time variable and append to bootstrap data
keep S_diff* _tscurve
append using data.dta
save data.dta, replace
restore
}

*load data of bootstrap survival curves
preserve
use data.dta, clear
keep S_diff* _tscurve
*drop cases with missing time variable
keep if _tscurve!=.
*reshape long S_diff1_ S_diff2_ S_diff3_ S_diff4_ S_diff5_, i(_tscurve) j(simulation)
*use percentile method to calculate upper and lower bound for each time point
forvalues i=1/5 {
bysort _tscurve: egen low`i'=pctile(S_diff`i'), p(2.5)
bysort _tscurve: egen high`i'=pctile(S_diff`i'), p(97.5)
}

*collapse to keep only one value for each time point and store data of results
collapse (min) low* (max) high*, by(_tscurve)
rename _tscurve _tscurve_boot
save data.dta, replace
restore

*merge bootstrap results with original data
merge using data
drop _merge


*the limited bootstrap sample gives a somewhat jagged confidence interval
*to be conservative, I use the outer bound for each time-period without a 
*change in survival curve difference
*this means that the confidence intervals will be smooth, although a bit to 
*large. The process is repeated for all five strata.
forvalues i=1/5 {
capture drop _change_id
gen _change_id=1 if Sdiff`i'[_n]!=Sdiff`i'[_n-1]
replace _change_id=sum(_change_id)
bysort _change_id: egen high`i'_n=max(high`i')
bysort _change_id: egen low`i'_n=min(low`i')
rename high`i'_n Sdiff`i'_uci
rename low`i'_n Sdiff`i'_lci
}
drop high* low*






*graph results for each stratum
twoway line Sdiff1 _tscurve if _tscurve<3650, sort c(J) lc(black) lp(solid) /*
*/	 || line Sdiff1_* _tscurve_boot if _tscurve_boot<3650, sort c(J J) /*
*/		lcol(black black) lp(shortdash shortdash) legend(off) xtitle("")/* 
		legend(order(1 2) size(small) region(lwidth(none)))  
*/		saving(order1, replace) title("1 previous crisis") yline(0, lc(gs6))
forvalues t=2/4 {
twoway line Sdiff`t' _tscurve if _tscurve<3650, sort c(J) lc(black) lp(solid) /*
*/	 || line Sdiff`t'_* _tscurve_boot if _tscurve_boot<3650, sort c(J J) /*
*/		lcol(black black) lp(shortdash shortdash) legend(off) xtitle("")/* 
*/		saving(order`t', replace) title("`t' previous crisis") yline(0, lc(gs6))
}
twoway line Sdiff5 _tscurve if _tscurve<3650, sort c(J) lc(black) lp(solid) /*
*/	 || line Sdiff5_* _tscurve_boot if _tscurve_boot<3650, sort c(J J) /*
*/		lcol(black black) lp(shortdash shortdash) legend(off) xtitle("")/* 
*/		saving(order5, replace) title("5 or more previous crisis") yline(0, lc(gs6))

graph combine order1.gph order2.gph order3.gph order4.gph order5.gph, ycommon b2("Days at peace", ring(0)) l2("Difference in survival probability")
graph export Beardsley_survival_diff_cox_atmean_strata.pdf, as(pdf) replace



